using System;
using System.Collections.Generic;
using System.Linq;

using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.LinearAlgebra;

namespace Atomic.Thermodynamics.StateEquations
{
	[Serializable]
	public class BirchMurnaghan : IStateEquation
	{
		private int n;
		private double e0, v0, b0, bp0;

		public BirchMurnaghan(Vector coefficients, double v0, int n = 2)
		{
			if (coefficients == null)
			{
				throw new ArgumentNullException();
			}

			if (n != 1 && n != 2)
			{
				throw new ArgumentOutOfRangeException();
			}

			Coefficients = coefficients;

			if (n == 2)
			{
				if (coefficients.Length == 4)
				{
					double bpp0;
					Transform(4, coefficients[0], coefficients[1], coefficients[2], coefficients[3], 0.0, 0.0, out e0, out v0, out b0, out bp0, out bpp0);
				}

				if (coefficients.Length == 5)
				{
					// Use equilibrium volume initial guess.
					double bpp0;
					Transform(5, coefficients[0], coefficients[1], coefficients[2], coefficients[3], coefficients[4], v0, out e0, out v0, out b0, out bp0, out bpp0);
				}
			}

			this.n = n;
			this.v0 = v0;
		}

		public static BirchMurnaghan Fit(IEnumerable<IVolumeEnergyPoint> points, int terms, int n = 2)
		{
			if (points == null)
			{
				throw new ArgumentNullException();
			}

			if (terms < 0)
			{
				throw new ArgumentOutOfRangeException();
			}

			if (n != 1 && n != 2)
			{
				throw new ArgumentOutOfRangeException();
			}

			if (terms == 4)
			{
				double e0, v0, b0, bp0, bpp0;
				Vector c = Solve4(points.ToArray(), out e0, out v0, out b0, out bp0, out bpp0);
				return new BirchMurnaghan(c, v0, 2);
			}

			if (terms == 5)
			{
				double e0, v0, b0, bp0, bpp0;
				Vector c = Solve5(points.ToArray(), out e0, out v0, out b0, out bp0, out bpp0);
				return new BirchMurnaghan(c, v0, 2);
			}

			throw new NotImplementedException();
		}

		public static Vector Solve4(IVolumeEnergyPoint[] points, out double e0, out double v0, out double b0, out double bp0, out double bpp0)
		{
			Vector x = SolveLeastSquare(points, 4);
			Transform(4, x[0], x[1], x[2], x[3], 0.0, 0.0, out e0, out v0, out b0, out bp0, out bpp0);

			return x;
		}

		public static Vector Solve5(IVolumeEnergyPoint[] points, out double e0, out double v0, out double b0, out double bp0, out double bpp0)
		{
			// Use v0 in BM4 as initial guess.
			Solve4(points, out e0, out v0, out b0, out bp0, out bpp0);

			Vector x = SolveLeastSquare(points, 5);
			Transform(5, x[0], x[1], x[2], x[3], x[4], v0, out e0, out v0, out b0, out bp0, out bpp0);
			return x;
		}

		private static Vector SolveLeastSquare(IVolumeEnergyPoint[] points, int terms, int n = 2)
		{
			int m = points.Length;

			double[] b = new double[m];
			double[,] a = new double[m, terms];
			for (int i = 0; i < m; i++)
			{
				b[i] = points[i].Energy;
				a[i, 0] = 1.0;
				for (int j = 1; j < terms; j++)
				{
					a[i, j] = Math.Pow(points[i].Volume, -n / 3.0 * j);
				}
			}

			// Least square fitting using pseudo-inversion.
			return SingularValueDecomposition.PseudoInverse(new Matrix(a)) * new Vector(b);
		}

		private static void Transform(int terms, double a, double b, double c, double d, double e, double v, out double e0, out double v0, out double b0, out double bp0, out double bpp0)
		{
			double b3 = b * b * b;
			double c2 = c * c;
			double c3 = c2 * c;

			if (terms == 4)
			{
				v = Math.Sqrt(-(4.0 * c3 - 9.0 * b * c * d + Math.Sqrt(c2 - 3.0 * b * d) * (4.0 * c2 - 3.0 * b * d)) / b3);
			}
			else
			{
				// Use supplied initial guess from BM4.
				v = RootFinder.Bisection(x => SolveVolume5(a, b, c, d, e, x), v * 0.9, v * 1.1);
			}

			double v2 = v * v;
			double v23 = Math.Pow(v, 2.0 / 3);
			double v43 = v23 * v23;
			double v83 = v43 * v43;
			double v113 = v83 * v;
			double v133 = v113 * v23;

			e0 = a + e / v83 + d / v2 + c / v43 + b / v23;

			b0 = 2.0 * (44.0 * e + 27.0 * d * v23 + 14.0 * c * v43 + 5.0 * b * v2) / (9.0 * v113);
			bp0 = (484.0 * e + 243.0 * d * v23 + 98.0 * c * v43 + 25.0 * b * v2) / (132.0 * e + 81.0 * d * v23 + 42.0 * c * v43 + 15.0 * b * v2);
			bpp0 = (4.0 * v133 * (27.0 * d * (22.0 * e + 7.0 * c * v43 + 10.0 * b * v2) + v23 * (990.0 * b * e * v23 + 7.0 * c * (176.0 * e + 5.0 * b * v2)))) / Math.Pow(44.0 * e + 27.0 * d * v23 + 14.0 * c * v43 + 5.0 * b * v2, 3.0);

			v0 = v;
		}

		private static double SolveVolume5(double a, double b, double c, double d, double e, double x)
		{
			double x3 = x * x * x;
			double x23 = Math.Pow(x, 2.0 / 3.0);
			double x53 = x * x23;
			double x73 = x53 * x23;
			double x113 = x73 * x23 * x23;
			return 8.0 * e / (3.0 * x113) + 2.0 * d / x3 + 4.0 * c / (3.0 * x73) + 2.0 * b / (3.0 * x53);
		}

		public double Energy(double v)
		{
			double x = Math.Pow(v, -n / 3.0);
			double e = 0.0;
			for (int i = Coefficients.Length - 1; i >= 0; i--)
			{
				e = Coefficients[i] + x * e;
			}

			return e;
		}

		public double Pressure(double volume)
		{
			throw new NotImplementedException();
		}

		public double BulkModulus(double volume)
		{
			throw new NotImplementedException();
		}

		public Vector Coefficients
		{
			get;
			private set;
		}

		public double EquilibriumEnergy
		{
			get
			{
				return e0;
			}
		}

		public double EquilibriumVolume
		{
			get
			{
				return v0;
			}
		}

		public double EquilibriumBulkModulus
		{
			get
			{
				return b0;
			}
		}

		public double EquilibriumBulkModulusPressureDerivative
		{
			get
			{
				return bp0;
			}
		}
	}
}
